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^ ■ ABSTRACT 

> ; 

^ ! A general method for constructing high order upwind schemes for multidimensional 

Q I magnetohydrodynamics (MHD), having as a main built-in condition the divergence- 

O ' free constraint V ■ B = for the magnetic field vector B, is proposed. The suggested 

^Z* • procedure is based on consistency arguments, by taking into account the specific op- 

Q^ ! erator structure of MHD equations with respect to the reference Euler equations of 

"^ I gas-dynamics. This approach leads in a natural way to a staggered representation of 

Ph' the B field numerical data where the divergence-free condition in the cell-averaged form, 

Q ■ corresponding to second order accurate numerical derivatives, is exactly fulfilled. To ex- 

^ ! tend this property to higher order schemes, we then give general prescriptions to satisfy 

^1 a (r + 1)*'^ order accurate V ■ B = relation for any numerical B field having a r*'^ 

l; ' order interpolation accuracy. Consistency arguments lead also to a proper formulation 

of the upwind procedures needed to integrate the induction equations, assuring the exact 
conservation in time of the divergence-free condition and the related continuity prop- 
el I erties for the B vector components. As an application, a third order code to simulate 
multidimensional MHD flows of astrophysical interest is developed using ENO-based re- 
construction algorithms. Several test problems to illustrate and validate the proposed 
approach are finally presented. 
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1. Introduction 

Many astrophysical plasmas, such as stellar (or 
galactic) atmospheres and winds, accretion disks and 
jets, can be described by the set of compressible mag- 
netohydrodynamic (MHD) equations with dissipative 
terms neglected, since kinetic effects of astrophysi- 
cal plasmas are quite small on dominant macroscopic 
scales. In these physical regimes, dynamical effects 
give rise to complex time dependent flows where lo- 
calized sharp modes like shocks and current sheets 
couple with distributed nonlinear waves. It is there- 
fore a main challenge to computational astrophysics 
to take properly into account both dynamical compo- 
nents. 

Centered finite differences or spectral schemes are 
well suited for smooth fields and can support discon- 
tinuities only by introducing enough viscous/resistive 
dissipation. In this way field discontinuities are rep- 
resented with a poor resolution and artificial heat- 
ing takes place. On the other hand, upwind schemes 
achieve shock-capturing and localized high resolution 
in a natural way. When discontinuous solutions are 
of main interest, second order (in time and space) 
schemes are usually adopted, since they reconcile res- 
olution with efficiency and stability needs. But in the 
general case, when coherent sharp field structures are 
embedded in a turbulent background, second order 
accuracy is no longer the optimal one, since the (im- 
plicit) numerical viscosity is still to high to resolve 
properly small scales motions. There are then com- 
pelling computational and physical reasons to develop 
higher order upwind schemes for MHD fiows. 

In recent years progress has been made in extend- 
ing Godunov-type schemes developed for the Euler 
system of gas-dynamics to MHD, with main empha- 
sis on the wave characteristic structure. In Brio & 
Wu (1988), and Roe & Balsara (1996) the problem of 
non strict hyperbolicity of the MHD system has been 
addressed by introducing proper regularity factors to 
renormalize the eigenvectors and assure their linear 
independence. The related problem of constructing 
the Roe linearized matrix for MHD case, has also been 
solved (Cargo & Galhce 1997, Balsara 1998a). Based 
on these achievements, second order upwind codes us- 
ing either Godunov's or Roe's method have then been 
constructed and tested, mainly for one-dimensional 
MHD problems (e.g. Ryu & Jones 1995, Zachary et 
al. 1994, Dai & Woodward 1994, Balsara 1998b). 

Specific new problems and limitations have to 



be considered in going to higher order and multi- 
dimensional MHD case. Upwind schemes are usu- 
ally constructed by first projecting fluid variables at 
each grid point on the space of characteristic vari- 
ables. The decomposition procedure allows to achieve 
a better resolution since interacting discontinuities 
of fluid variables become uncoupled in the space of 
characteristic variables. This technique is usually 
adopted also in existing MHD codes, but there is 
no clear evidence that it can work even for higher 
order schemes (Barmin et al. 1996). Moreover, the 
computational cost to project field variables onto the 
MHD seven-component characteristic space may be- 
come prohibitive when moving to higher order and 
higher dimensional schemes. Therefore, as already ex- 
perienced in the context of numerical gas-dynamics, a 
search for high order shock-capturing schemes where 
no characteristic decomposition is needed and where 
attention is shifted more to a vanishing viscosity 
entropy satisfying model equations, rather than on 
approximate Riemann solvers, appears to be more 
promising. 

A second important issue in multidimensional MHD 
schemes comes from the need to satisfy the divergence- 
free condition of the magnetic field vector. This prop- 
erty is a crucial one for two main reasons (Balsara & 
Spicer 1999): 

• the conservation form of MHD equations for 
energy and momenta is based on the implicit 
V • B = condition, and 

• all the topological aspects of magnetic field lines 
which are relevant to critical MHD phenomena, 
like reconnection, heavily rely on this condition. 

On the other hand, this specific property has no 
easy representation in a numerical framework, like 1- 
D Godunov-type schemes, designed to handle com- 
pressive modes and shocks. This longstanding prob- 
lem has been addressed by many authors and several 
recipes have been proposed and experimented so far. 
Depending on the adopted methodology, these works 
can be broadly classified into three main categories: 

1. Many MHD codes are constructed by simply 
extending to higher dimensions 1-D Riemann 
solvers using a directional splitting technique, 
as for the Euler system (Zachary et al. 1994, 
Ryu et al. 1995, Balsara 1998b, for second or- 
der Godunov-type schemes; Jiang & Wu 1999, 



for S*'' order WENO scheme). In this approach 
the V • B = condition breaks down, of course, 
and some correction step has then to be apphed. 
Following Brackbill & Barnes 1980, a cleaning 
procedure is usually carried out by solving a 
Poisson equation, which is equivalent to add 
a new (elliptic) equation to the original hyper- 
bolic MHD system. As an empirical recipe, this 
method is by no means optimal and may lead 
to inconsistencies. In particular, the numer- 
ical derivatives appearing in Poisson equation 
have no clear relation with the upwind deriva- 
tives of the base MHD system and the boundary 
conditions become indeterminate for nontrivial 
boundary- value problems. 

2. In the Powell (1994) approach it is first pointed 
out the formal difhculty of applying 1-D Rie- 
mann solvers to the multidimensional case, since 
the 1-D MHD mode eigenspace, having seven 
components, is not of full rank for the 2-D case, 
where an eight-component state vector is in- 
volved. Therefore, variations of magnetic field 
components appearing in the V-B operator can- 
not be represented by Riemann solvers based on 
the 1-D eigenspace. This undoubtedly correct 
premise led the author to propose a modifica- 
tion of the MHD equations by adding a "V • B 
mode" , propagating with the local flow speed. 
In this way the hyperbolic character of the MHD 
system is surely retained, at the price of sup- 
pressing the divergence-free property. This ap- 
proach appears to be highly questionable, of 
course, since important physical properties of 
the MHD equations, and especially magnetic 
field topologies, are clearly lost. 

3. In the present work we have taken as a main 
starting point all those references attempting 
to design upwind schemes where a numerical 
divergence-free condition works as a build-in 
property (Evans & Hawley 1988, DeVore 1991, 
Stone & Norman 1992, Dai & Woodward 1998, 
Ryu et al. 1998, Balsara & Spicer 1999, among 
others). In all these works, the introduction 
of the magnetic vector potential or the equiv- 
alent conservative formulation of Stoke's theo- 
rem lead to represent magnetic field components 
at staggered collocation points, and a numerical 
V • B — relation follows as an algebraic iden- 
tity. Moreover, when induction equations are 



properly formulated in terms of the staggered 
fields, conservation in time of the divergence- 
free property is also assured. 

In the cited works, however, some main ques- 
tions are still left open. These are essentially 
related to a persisting duality between staggered 
magnetic components evolving in the induction 
equations and the same components, now collo- 
cated at node points (or cell centers) as other 
fluid variables, entering the Riemann solver pro- 
cedures. Several recipes based on interpolation 
have been suggested to relate cell centered and 
staggered fields. However, as widely discussed 
in the Dai & Woodward (1998) paper, the cell 
centered field components do not preserve, in 
general, the original divergence-free property, 
unphysical magnetic monopoles still arise and 
their sizes seem to depend on the adopted in- 
terpolation schemes. 

A related question concerns how upwind fluxes 
in the induction equations have to be formu- 
lated, since 1-D Riemann solvers for density, 
momentum and energy equations have no straight- 
forward extension to them, when staggering is 
adopted. Again, many different empirical solu- 
tions have been proposed, which hardly can be 
compared and evaluated as long as only quali- 
tative numerical tests are at disposal. 

In the present paper we propose some general an- 
swers to these questions, by showing that consistency 
arguments are sufficient to envisage the main rules 
to adapt upwind schemes designed for Euler equa- 
tions to the MHD case. Consistency requires that 
the specific operator structure of the MHD system 
and the related magnetic field properties have to be 
preserved by discretized equations and upwind proce- 
dures. In this way, different schemes and their high 
order extensions can be designed, all assuring a nu- 
merical divergence-free condition as well as the re- 
lated uniqueness and regularity of magnetic field lines. 
On the same ground, existing MHD codes and pub- 
lished numerical results can be evaluated on a more 
appropriate framework. 

The plan of the paper is as follows. In Sect. 2 
the general formulation to discrctize Euler and MHD 
equations in conservation form is reviewed, with em- 
phasis on the differences in space operator struc- 
tures and on the related numerical representation. 
In Sect. 3, the kinematical properties of discontinu- 



ous, divergence- free magnetic field are first analyzed 
to represent numerical data and then used to con- 
struct appropriate upwind flux formulas (or approxi- 
mate Riemann solvers) for the MHD equations. The 
proposed formulation is then also compared to re- 
cently published schemes. In Sect. 4 a code for 2-D 
systems, based on third order ENO-type reconstruc- 
tion procedures and on the simple Lax-Friedrichs flux 
upwinding, is presented. Sect. 5 is devoted to numer- 
ical test problems to add confidence and validation 
of the proposed approach and conclusive remarks are 
briefly given in Sect. 6. 

2. Euler versus MHD systems 

To underline differences between the Euler and 
MHD systems relevant to numerical discretization, we 
first briefly review some of the main points character- 
izing upwind schemes for the Euler equations. As a 
general framework, we consider here the flux vector 
splitting (FVS) formalism (van Leer 1982, Chen & 
Lefloch 1995) and the high order reconstruction tech- 
niques based on polynomials (Shu 1997). For ease of 
presentation we treat only the 2-D case in cartesian 
geometry, being both the 3-D case and curvilinear ge- 
ometries just straightforward extensions. 

2.1. Upwind schemes for Euler equations 

The equations of gas-dynamics in two spatial di- 
mensions constitute a system of to = 5 conservation 
laws: 

dtU + dSin)]+dy[giu)]^0, (1) 

where u — [p, q, e]'^ is the state vector of conservative 
variables, and 

f = [qx^v^qx + p,Vxqy,Vxq^,Vx{e + p)f , 

g = [qy,Vyqx,Vyqy + p,Vyq^,Vy{e + p)]'^ , 

are the corresponding flux vector functions. Here p is 
the mass density, q = pv the momentum associated 
with the flow velocity v, e the total energy per unit 
volume and p = (7— l)[e — q- v/2] the gas pressure 
for a 7— law equation of state. 

A basic property of the Euler system is that each 
Jacobian matrix, A2,(u) = 9uf(u) and Aj,(u) = 
9ug(u), has a set of m real eigenvalues {A''(u)} (s = 
1, 2, . . . , m), and a corresponding complete set of right 
{Rs(u)} and left {Rj^(u)} eigenvectors, at every 
point u (hyperbolicity properties). Physically rele- 
vant solutions to system (nh are selected by imposing 



the admissibility condition 

^t[pF{s)]+V■[pF{s)^^]<0, 



(2) 



where F{s) is any smooth function of the specific en- 
tropy s{p,p) (Harten et al. 1998). 

Numerical schemes for system (|^) use the following 
consistency conditions as general guidelines (Tadmor 
1988): 

• the conservation form, assuring that numerical 
solutions capture correctly weak solutions; 

• the entropy inequality, to be preserved by the 
discretized entropy functions. Upwind schemes 
are then designed to have (implicit) numerical 
viscosity compatible with relation (g). 

In the scnii-discrete formulation, appropriate for higher 
r > 2 order schemes, space operators are approx- 
imated (for a fixed time t) on a Nx x Ny dimen- 
sional grid with node points Nj_k = {xj,yk), where 
Xj = hxj (3 = 0, 1, . . . , TVa; - 1) and yu = hyk [k = 
0, 1, . . . , iV^y — 1); here hx and hy are the constant grid 
sizes along each direction. The point values formula- 
tion based on {u^-jt} data leads then to the conserva- 
tive scheme 






dt 



Ij-l/2,fcJ 



1 



lj.k+l/2 



?j",fe-l/2j7 



(3) 



where ij+i/2,k and g,j^k+i/2 denote the numerical flux- 
vector functions needed to approximate the corre- 
sponding flux derivatives to a given order r. A numer- 
ical approximation is then characterized essentially by 
the way fj+i/2.fc f^nd g,j^k+i/2 are evaluated for a given 
set of {uj.fe(i)} data. Time integration can then be 
performed by appropriate Runge-Kutta or equivalent 
stable discretization schemes (Shu & Osher 1988). 

Modern higher order shock-capturing schemes gen- 
eralize first order Godunov scheme by following two 
main steps (Harten et al. 1987): 

1 . a reconstruction phase to recover variable values 
at grid points where flux derivatives have to be 
computed; 

2. an upwind phase, where the Godunov method 
for a scalar conservation law in one dimension 
is extended to the m > 1 components system in 
higher dimensions. 



As far as item 1 is concerned, we remind here some 
basic points relevant to the following sections and to 
the actual code structure, to be presented in Sect. 4. 

Any one-dimensional piecewise smooth function 
w{x), defined by cell averaged data 



1 



'w{x)dx, 



can be approximated by uniform (r — 1)*'' order piece- 
wise polynomials Pj{x;w), which must have the con- 



servative property Pj 



Likewise, a w(x) func- 



tion defined by grid-point data {wj = w{xj)} can be 
approximated by interpolation polynomials Pj{x]w) 
defined as Pj{xj;w) — Wj. Here we denote with 
R[x; W] (or R[x; w]) the corresponding polynomials set 
{Pj} to reconstruct w{x) at any x point. 

At points of discontinuity, R[x; ■] has to satisfy 
definite non-oscillatory constraints, for accuracy and 
stability purposes. Standard references are provided 
by linear polynomials based on minmod limiters to 
preserve monotonicity of data (MUSCL scheme: van 
Leer 1979; TVD scheme: Harten 1983) or by higher 
r > 2 order polynomials based on ENO procedures 
(Harten et al. 1987, Shu & Osher 1989) having weaker 
(essentially non-oscillatory) monotonicity properties. 
Piecewise polynomials approximate the w{x) function 
at any cell boundary point 2:^+1/2 by a two-point left- 
right {w^-'"\ w^-^^) value, where 



[^^''^], + l/2 



P. 



j(Xj + y2)=W[X-^^^^) 



0{K) 



\w 



(RY 



Xj < X < Xj+i range, with 



j+1/2 = Pj+i{xj+i/2) = ^(x+^^/J + Oihl), 

for smooth functions. For w{x) having a discontinu- 
ous fc*^ derivative in the 
k < r, the accuracy order becomes 0{h^'^^) 

The reconstruction procedures can be extended to 
the 2-D functions in Eq. (P by assuming that the 
scalar components u{x,y) (and hence f[\i{x,y)] and 
g[\i{x,y)\) are piecewise smooth along each coordi- 
nate. In this way, calling Cj^k the 2-D cell cen- 
tered at the node point Nj^k, the functions u(x,y) 
can be reconstructed at the cell edges (xj+i/2,j/fc) 
and (a^j, yfc+1/2) by using respectively the 1-D op- 
erators R[x;uk] and R[y]Uj\. The same functions 
may also be reconstructed at a cell corner Pj^ = 
{xj^i/2, yk+1/2) by using the 2-D compound operators 
R[xj+i/2]R[yk+i/2]'^]]- Actually, space discretization 
in Eq. dl^) involves only 1-D reconstructions, one for 
each direction. In fact, fj+i/2,k is defined at the point 



X = Xj+if2, for fixed y — yk, where the argument vari- 
ables have the two-state a;— wise reconstructed values 



\u 



(L, 



, u^^'']k', likewise, gj^k+1/2 is defined at the point 



y = yk+1/2, 



for fixed x 



where the argument vari- 



ables have the two-state y— wise reconstructed values 

Let us now turn our attention to item 2. In the 
FVS formalism, to represent a fiux variation, say in 
the X coordinate, the vector f (u) is decomposed as 

f(u) = i[f(+) + f(-)], f(±) - f(u) ± D,(u) • u, 

for states u around a given reference (constant) state 
u. The fiux vectors f'^^-' have Jacobian matrices 
Ax (u) ± T>x (u) with only positive/negative eigenval- 
ues. The matrix Da; is then required to be real (sym- 
metrizable) and positive, 'Dx{ii) > |Aa:(u)|, where 



|A.(u)|=5][R 



= V[R,|A,|R;i]i 



For given Ui — u*^-^'\u2 — u^^--") reconstructed values 
at the point 2:^+1/2 and for fixed k index, flux splitting 
allows to define the numerical flux 



1 



f(Ul,U2)j + i/2,fe = 
[f(ui) + f(u2) - Da;(u) • (U2 - Ui)]j-+i/2,fc, 



(4) 



which has upwind properties, namely it is a two-point 
vector function non-increasing in the first argument 
and non-decreasing in the second argument. The ref- 
erence state u is given by the Roe average or by the 
simpler u — (ui + U2)/2 arithmetic average, in a way 
to assure consistency f (u, u) = f (u) and continuous 
dependence on data. 

Standard references for the scheme in Eq. (H) are 
given by either the approximate Riemann solvers of 
Roe-type or of Godunov-type , where the correspond- 
ing matrix Dj, has the formal property 

D,(u) = |A,(u)|+0(|u(^-)-u(^-)|) 

entailing a minimum of numerical viscosity compati- 
ble with the entropy law. On the other hand, a maxi- 
mum of numerical viscosity is achieved by the (global) 
Lax-Friedrichs (LF) flux, where the D^, matrix re- 
duces to the simple diagonal form 

D2;(u) = q;I, a — maxii[maXs\Xsii^)\], 

and in which Riemann characteristic informations are 
averaged out. Note that the flux formula Eq. (||) 



can be interpreted either as an approximate Riemann 
solver based on local linearization, or as a discrete 
approximation of the associated viscosity model equa- 
tion. In fact, a first order approximation of the flux 
splitting is the discretized representation of the vis- 
cous flux 



fi;(u) = f(u) - hx'Dx ■ dxU, 



(5) 



which provides a link between the entropy condition 
and the numerical viscosity associated to the dissipa- 
tion matrix hj-Ty^. 

In 2-D problems, the semi-discrete formulation of 
Eq. (0) and the independence of the Jacobian ma- 
trices (Aa;, Ay), allow one to represent the numerical 
fluxes f,+i/2,fe (for fixed yk) and gj,k+i/2 (for fixed Xj) 
through independent upwind procedures, constructed 
with the matrices T>x and Dy, respectively. In this 
way, the fiux formula for g is given by 

g(ui,U2)j,fe+i/2 = 
-[g(Ui) -I- g(u2) - Dy(u) • (U2 - Ui)]^- fc+i/2, (6) 

where Ui = u*^'^) U2 = u^'-'^', are the reconstructed 
values at y = ?/fe+i/2j for fixed x = Xj. The cor- 
responding viscosity form of the numerical fiux g is 
given by 



Svi^) = g(u) - hyBy ■ dyU. 



(7) 



The construction of upwind fluxes in Eqs. (0) and 
(|^), based on the A^^ matrix eigenspace at the point 
(a;j_|_i/2,2/fe) and on the Aj, matrix eigenspace at the 
point {xj,yk+i/2): is usually referred to as a direc- 
tional splitting setting. This procedure is consistent 
with the divergence form of space operators in Eq. (|l|) 
and implies that the space derivatives are obtained by 
summing the two flux differences computed both at 
the same time t. In the Strang-type formalism, which 
is widely adopted in the numerical astrophysics com- 
munity, a splitting procedure is also applied to the 
time evolution operators, by constructing the updated 
solution of a 2-D problem as a sequel of independent 
1-D problems, one for each direction, in turn. In 
the Euler system directional splitting or time split- 
ting procedures give (formally) equivalent results, at 
least for second order schemes. In the MHD case, 
however, this formal equivalence is definitely lost, as 
it will be discussed in the following. 



2.2. The two-dimensional MHD system 

The set of MHD conservation laws cannot be con- 
sidered as a simple extension of the Euler system, 
with just a higher number of state variables. Ac- 
tually, while the conservation form, the entropy law 
and the general hyperbolic properties are maintained, 
some specific differences related to the structure of the 
space operators have to be considered. 

In fact, the MHD system can be viewed as com- 
posed by two coupled subsystems, the first one con- 
taining space operators in the divergence form as in 
Eq. (n^) evolving density energy and momenta, and 
the second one, specific to the magnetic field evolu- 
tion, containing space operators in the curl form. In 
both subsystems the V • B = property of the vec- 
tor magnetic field enters in a substantial way and has 
then to be considered as a new consistency condition 
for numerical discretization. 

By specializing again to 2-D systems, the MHD 
equations are given by 



ftu + a,[f(w)]-f9j,[g(w)]=0. 



(8) 



for the six-component state vector u = [p, q, e, B^]'^ , 
coupled with the induction equations 

- dtBx + dyiliw) = 0, dtBy + dMi^) = 0, (9) 

for the (poloidal) vector field B = [Bx,By]'^. We 
denote the overall eight-component state vector as 
w = [u,B]"^. The flux vectors in Eq. (B) are given 

by 

f (w) = [Qx, Fx^x, Fx,y, Fx^z, Ex, Gx,z] , 
g^Wj ~ [%, -Ty.xT -t'y.y, -t'y,ZT Ey,Uy^z\ j 

where, for indexes i,j = x,y, z: 

F^j = ViQ-i - BiBj + H^jj, 

E,=v^{e + Ii)-B,{vjB,), 

Gij = ViBj -VjBi, 

in which the relations Fi_j = Fj^i and G^j = —Gj,i 
clearly hold. Here H — p + {BiBi)/2 and p = (7 — 
l)[e — {qiVi)/2 — {BiBi)/2] are respectively the total 
and the gas pressures. The common flux function of 
Eqs. (||) is defined as J7 = Gx,y = —Gy^x = VxBy — 

VyBx- 

For given B field, the subsystem (||) as an Eu- 
ler form, with independent Jacobian matrices of full 
TO = 6 rank, so that the upwinding procedures based 



on directional splitting of the previous section can 
be extended. Differences arise, however, for the in- 
duction equations. In fact, subsystem (^ is gener- 
ated by a unique flux function and has then only a 
one-dimensional eigenspace. This formal property is 
clearly related to the V • B = condition, as can 
be better evidenced by introducing a vector potential 
representation of the (poloidal) B components (here 
A = Az) 

B, = dyA, By = -d^A, (10) 

For given smooth {Bx,By) fields, A{x,y) always ex- 
ists as a one-valued differentiable function. For dis- 
continuous fields, Eq. ( p^ still holds in weak form, 
implying that A(x, y) is at least (Lipschitz) contin- 
uous along each coordinate. On the other hand, for 
given A{x,y,t), the system given by Eqs. (S) is fully 
equivalent to the one-component evolution equation 

dtA - r!(w) = 0, (11) 

coupled with Eqs. (p^. 

In a formal setting, for given state variables u, 
Eq. (O) is a Hamilton- Jacobi equation (Jin & Xin 
1998), and Eqs. (O) represent the associated hyper- 
bolic system. An important property of Eqs. (O) 
is that A{x,y,t) is continuous at all times and only 
discontinuities in its first derivatives may develop. 
Therefore, field lines defined by the isocontours A{x, y) 
const are allowed to have corners, but not jumps. 

The overall MHD system can then be viewed as a 
coupled system of a Hamilton- Jacobi equation and of 
a set of conservation laws in the Euler form. One con- 
sequence is that the Jacobian matrix A^; correspond- 
ing to the [f , VlY' vector flux is only of to = 7 rank and 
can represent characteristic modes of variables w^. = 
\vL{x,),By{x,)Y' , while the independent matrix Aj,, 
corresponding to the [g, — f2]"^ vector flux, can repre- 
sent variables Wy = [u(,2/), i33:(,?/)]^. It is evident 
that the missing degrees of freedom [Bx{x, ), By{, y)] 
are not evolutionary and cannot have a characteristic- 
based representation. 

A numerical schemes preserving these general prop- 
erties must then be be characterized by the following 
points: 

1. A numerical V • B = condition and its conser- 
vation in time imply that the induction equa- 
tions (P) have to be discretized using a unique 
flux function ri("w) located at common points. 
This entails necessarily a staggered collocation 
of the magnetic field scalar components. 



2. A divergence- free magnetic field is fully equiv- 
alent to its representation via a numerical vec- 
tor potential and likewise the evolution equa- 
tions (0), discretized as in item 1, can always 
be integrated via the scalar Eq. (fl]). 

3. Relevant to the reconstruction and upwind steps 
is that the magnetic field components are at 
least continuous along the respective longitu- 
dinal coordinates, while discontinuities, to be 
related to the MHD characteristic modes, can 
occur only along the respective orthogonal co- 
ordinates (see below). 

3. V • B = preserving upwind schemes for 
MHD equations 

In this section we concentrate on the correct collo- 
cation and reconstruction step for the numerical mag- 
netic field data, and then on the upwind flux formula- 
tion for the induction equations, in order to preserve 
the peculiar features of the MHD system as outlined 
just above. 

3.1. The reconstruction step 

While for given data {ui j} of u variables in Eq. (ph 
the reconstruction procedures follow the same lines as 
in the Euler system (0) , for the field B it is necessary 
to take into account the V • B — condition as a new 
kinematical constraint. For general piecewise smooth 
fields this condition is expressed in integral form by 



Bndl = 



dc 



(12) 



for any cell C, where i?„ = B • n and n is the unit vec- 
tor normal to the boundary line dC . By first choosing 
a cartesian cell with sides [2e, hy], the following conti- 
nuity condition for the y— averaged Bx{x) component 
at any point x comes out: 

B,{x + e)-B^{x-e)^0{e). 

The same argument leads to the continuity of the 
X— averaged By{y) field at any point y. 

If Eq. ( p^ ) is then integrated over a computational 
cell Cj^fe, one has 



hx[By{yk+i/2) - By{yk-i/2)]j = 



(13) 



where [Bx]k is the y average on the vertical cell side 
centered on yk and [By]j the corresponding x average 
over the horizontal cell side centered on Xj . 

Continuity conditions and Eq. (O) are thus the 
main ingredients to construct at any point a divergence- 
free numerical field. In particular, the continuity 
condition B^ix^) = Bx{x~) allows one to locate 
the [Bx{x)]k field at cell boundary points {a;j_|_i/2}, 
where all the other variables [wx{xj^i/2)]k are rep- 
resented by reconstructed two-state (left-right) val- 
ues. This property can be expressed in a formal 

way by setting [B^ ' ]j+i/2,k = [B^ " ]j+i/2,k- Cor- 
respondingly, By{y) can be located at cell boundary 
points {yk+1/2} and point values can be interpreted 

as i^y \k+i/2 - i^y Jfc+1/2- 

The reconstruction step for the {Bx,By) fields 
along the respective transverse coordinates leads to 
the point values 

[Bx{y)]j+i/2 = R[y;Bx], [By{x)]k+i/2 = R[x;By], 

which have relevance for upwind computations. In 

fact, at the yk+1/2 Point, J5a;(2/fe+i/2) = [Sp^si''^^] 
is a two-state variable, and can have only a transverse 
discontinuity line with a SyBx = [Bx — Bx ] jump. 
Similar arguments apply to the reconstructed values 
By{xj+i/2) = [By ,By ] at the x = Xj+1/2 bound- 
ary point, allowing only for a transverse discontinuity 

As already noticed by Evans & Hawley (1988), a 
main property related to condition (O) is that the nu- 
merical data {Bx, By} allow one to construct a unique 
(continuous) numerical potential A{xj^i/2, yk+1/2): 
located at the cell corners Pj.k- It is evident that also 
the reverse condition holds true, by first introducing a 
numerical continuous function A{xj^i/2, yk+1/2) and 
then defining the averaged magnetic field components 
through Eqs. (0): 

[BxiXj + l/2)]k = —[AyA{Xj + i/2)]k, 



[Byiyk+i/2)h = -i-[^ccA{yk+i/2)h, (14) 

i^x 

where (A^., Aj,) denote the usual (undivided) centered 
finite differences on the first and second coordinate in- 
dex, respectively. In this way the divergence-free con- 
dition, Eq. (p^, is identically satisfied by the commu- 
tativity of the Ax and Ay linear operators, whereas 
the continuity conditions follow from the continuity 



of A{x, y). These arguments show, in particular, that 
the staggered collocation for {Bx,By} data is by no 
means a numerical trick but arises in a consistent way 



from Eq. (|l3|) or Eqs. (|lj). 

Equation (O) gives the cell average of the V-B = 
condition, and thus it is an exact law for second order 
accurate schemes since Bx = Bx + 0{hV) at any point 
{xj+i/2,yk) and By = By + OQi^.) at the correspond- 
ing staggered point. Now, if higher order approxima- 
tion of point- valued Bx and By fields were recovered 
using independent reconstruction steps based on Bx 
and By data, a numerical V • B of any size could arise, 
in general, since 1-D reconstruction operators do not 
commute. To overcome this main difficulty, a differ- 
ent strategy has to be adopted, by first reconstructing 
accurate first derivatives based on the vector poten- 
tial representation and having a numerical V • B = 
relation as a build-in property. As an illustration, in 
the following we consider third order interpolations, 
but extensions to higher order can be easily pursued. 

We first notice that cell averaged Wj and point val- 
ues Wj data of a given w{x) function are related by 

0{h^), where 71 = 1/24 



,(2) 



^iV'i\w)j 



w 

and "Dx"^' denotes a nonoscillatory numerical second 

derivative along the x coordinate. To the same ac- 

(2) 
curacy, the inverse relation Wj — Wj — jiVx {w)j 

approximates point values using averaged data. This 

algorithm, now applied to the values Wj+1/2 at cell 

(2) 

interfaces, w)j+i/2 = [w - 71 Pi (w)]j+i/2, gives the 
numerical 'w{x) primitive function whose two-point 
difference [Axw]j/hx constitutes a third order ap- 
proximation of the dxw{x) first derivative at Xj. Let 
then apply this recostruction step to approximate the 
primitive A of the magnetic potential A{xj+i/2, yk+1/2) 
In the 2-D (a;, y) plane we have 

[A]o + l/2.k+l/2 = [A- 7l(2?i'^ + I?.t'^)^]j + l/2,fe+l/2, 

and the numerical magnetic field components are 

1 
h 



[Bx\]+i/2.k - T-l^yA]]- 



1/2, fc, 



[By\j^k+l/2 — ~T-[^xA]j^k+l/2- 



(15) 



By definition, the difference [AxBx]/hx gives a third 
order approximation of the dxBx first derivative at 
the node {xj, yk) point, and [AyBy\/hy gives the cor- 
responding approximation of the dyBy derivative with 
the same accuracy and at the same point. We no- 
tice that no left or right derivatives are defined along 



the longitudinal coordinates, thus these numerical ap- 
proximations are unique. Moreover, one easily verifies 
that V • B = 0, now in the point-valued form, is ex- 
actly fulfilled due to the commutativity of the A^Ay 
operator. The key point here is that to higher orders 
only the primitives [Bx{x,), By{,y)] can be recon- 
structed directly using a common magnetic potential 
function A(x, y), but notthe (B^, By) functions them- 
selves (those entering the fluxes, where divergence- 
free fields are actually needed). To achieve this, one 
needs a further computational step, that is to solve 
the (now implicit) relations 



[B,-jiV^^^B.. 



x\j+l/2M 



[B: 



x\j + l/2,k: 



[By - jiV(y^^By]^^k+i/2 = [By]j,k+i/2, (16) 



where on the left hand sides appear the primitives, 
and hence the derivatives, defined in terms of the 
(unknown) field point values, while on the right hand 



sides are the source terms {Bx,By), given by Eqs. ([151). 
The numerical fields (B^^By) defined by these equa- 
tions and the corresponding (longitudinal) deriva- 
tives are third order approximations but satisfy the 
divergence-free condition exactly. In practice one can 
solve Eqs. (|l^) by some explicit iterative algorithm, 
since each operator (1 — 7il>(^)) is clearly invert- 
ible and a fourth order accurate, at least, V • B = 
0{h'^,hy) condition can then be easily satisfied (see 
Sect. 4). 

This completes the main proof for the reconstruc- 
tion step, needed to represent magnetic field point 
values in the momentum and energy equations, where 
longitudinal derivatives and hence a V • B = condi- 
tion has to be satisfied to avoid numerical monopoles. 
Moreover, for given [Ba:]j+i/2,k and [By]j ,,+1/2 data 
it is possible to get interpolated values at other col- 
location points, where these field components act as 
independent variables and no divergence-free condi- 
tion is then required. 

3.2. The upvirind procedures 

Let us now consider the first set of MHD equa- 
tions (H) discretized as in Eq. (g): 

duj^kjt) _ 
dt ~ 

~ T~[^j+l/2,k ~ ^j-l/2,k\ ~ T~[Sj,k+l/2 ~ gj,A;-l/2]j 

fix ify 

(17) 



where now the flux functions [f (w), g(vv)] depend on 
the eight-component vector w = [u,B]. The upwind 
flux based on the A^, characteristic eigenspace has the 
form (consult Eq. (0)): 

[f(w„i?,)],+i/2,fc = i[f(wi^'),i3,) +f(wi^'),B,) 

-Dr6)(w).(wi«-)-wi^^))],-+i/2,fc (18) 

where now upwind properties involve only the vari- 
ables Wa; — [\i^By\. In the same way, the upwind 
flux based on the Ay characteristic eigenspace has 
the form (consult Eq. (H)): 

[g(w„i3,)],-fe+i/2 = i[g(w(^^),S,) +g(w('^),S,) 

^D(i-6)(w).(w(^^)-w('^))],-,fe+i/2, (19) 

where this time the upwind properties involve only 
the variables Wj, = [u, B^]- 

To second order approximation, the numerical flux 
needed to compute space derivatives in Eq. (|l^) are 
given by the flux values of Eqs. ([l8| , [l9| ), whose argu- 
ments yvx and Wy are second order interpolated vari- 
ables. In particular B^ = B^) in the ^{'WxtBx) flux 
and By = By in the g(wj,. By) flux. 

On the other hand, in classical second order schemes 
for the Euler equations the numerical flux f is ex- 
pressed using cell centered variables as 

fj-Hl/2,fc = 1^[^{^3 + I,k) + f(Wj,fc)] 

-i[Dli-«)(w).(w(^^)-w^))Ui/,,, 

and in a similar way for the g flux. However, this flux 
representation cannot be extended to the MHD case 
since cell centered {B^, By) flelds are not related by a 
divergence-free condition. The same remark applies 
to higher r > 2 order schemes where the numerical 
flux reconstruction is based on cell centered values 
(Shu & Osher 1989). 

A second consequence of MHD structure in the 
evolution equations (^ is that the u state vector 
has to be integrated in time by summing flux deriva- 
tives evaluated at the sam,e tim,e t, when the implicit 
V • B = condition holds. In those schemes where 
time integration is performed by a Strang-type split- 
ting procedure, flux derivatives and hence terms con- 
taining dxBx and dyBy are necessarily summed at 



different time steps and the required dxB^ + dyBy 
cancellation never occurs. 

To summarize, higher order upwind schemes de- 
veloped for Euler equations can be extended to the 
MHD sub-system ([l7|), with the proviso 

1. flux derivatives have to be computed using flux 
values and hence B field data directly collocated 
at staggered (i.e. cell boundary centered) points 
and not at the cell centered points and, 

2. the same derivatives along the two directions 
have to be computed at the same time, thus 
avoiding time-splitting techniques. 

As already anticipated, the second set of MHD 
equations, given by the induction equations (||) for 
the magnetic poloidal field components, needs a par- 
ticular treatment. In the proper, divergence-free pre- 
serving discretized form, these equations are given by 

j^[By{t)]^^k+i/2 = -^[A^niwp)], 

^[B.(i)],+i/2.fc = ^[Ayn{wp)], (20) 

where now w = [vx,Vy,Bx,By]'^ denotes only the 
variables which are arguments of f2. In a fully equiv- 
alent form using the vector potential representation, 
one has 



Di7)(w).5,w,](y),+i 



/2, 



(22) 



dt 



[A{t)]p - f}(wp). 



(21) 



to be coupled with Eqs. (|lj). In both formulations a 
common flux function ri(w), located at the cell corner 
point P = [a;j+i/2, J/fc+i/2]i has to be evaluated. 

In order to single out a consistent numerical flux 
function in Eq. (|2l| ) , one has to take into account that 
n(wp) is now a four-state function 

and upwind rules involve necessarily both A^ (w) and 
Ay{w) matrix eigenspaces evaluated at a common ref- 
erence state [w]p. To construct a proper 2-D Rie- 
mann flux formula, we first consider the two limiting 
cases in which a propagating discontinuity may be 
described by a 1-D Riemann flux formula. 

1. For a discontinuity front perpendicular to the x 
direction, where w^'^^ ~ w^'^' , one has 

[f^.(y)Ui/2 = ^[f^(w(^')) + f^(w(^-)) 



»(7) 



where Dx denotes the seventh row vector com- 
ponent of T)x matrix and Sx'Wx = (wi; — 
wi^'^). The flux in Eq. (JH) is extended to 
the range yk l£ y l£ 2/fe+i where the vari- 
ables w(°')(y), a = L,R are continuous. We 
remind here that the upwind formula ( p2[ ) can 
be derived from the flux splitting introduced in 
Sect. 2.1, in the form 

n'^^^'i = n{w)±B'J\^)-wx, (23) 

which is equivalent to a local linearization of the 



n{'w{x)) flux around the x 



^1/2 point. 



2. For a discontinuity front perpendicular to the y 
direction, where now w^^'^ = w^^'^ the approx- 
imate Riemann solver reduces to 



[f^.(^)Wl/2 = ^[f^(w('^)) + f^(w(^^)) 



■Dl,^^(w)-(5j,Wy](a;)fc+i 



1/2, 



(24) 



^m 



where Dy denotes the seventh row vector com- 
ponent of the Dj, matrix and SyWy — {wy — 

Wy '). The X range involved is now Xj < x < 
Xj+i where the variables w'^'''^{x), b = L,R are 



continuous. Again, the upwind formula (24) 
may be derived from the splitting 



n 



(.±) 



f}(w) T D^'^w) . w. 



(25) 



which represents a local linearization around the 
y = 2/fc+i/2 point. 

In the general case where a discontinuity front 
crosses the computational cell centered on P, an ap- 
proximate Riemann solver can be obtained by intro- 
ducing a 2-D flux splitting. For that purpose, we 
decompose each fl'^^''^ components in Eq. (23) along 
the y direction, with the requirement to have the 
same form of the symmetric decomposition of the 
ri'^^'^'^ components in Eq. ( P3| ) along the x direction. 
This compound flux splitting, when interpreted as an 
approximate Riemann solver with local linearization 
(thus implying to neglect 0{Sx'WxSyWy) terms), re- 
sults in the four-state flux formula 

n{wp) = 

i[f}(w(«'«)) + r!(w(^'^)) + f^(w(^-«)) + f^(w(^'^))]p 
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hBP{^)-S^W,~'D^J\^)-SyWy]l 



(26) 



where (a, b ^ R,L) 



w^-^^) - 



y 2 ^ 



.("^^i)^ 



,(i,b) 



), 



andw= (w(^-")+w('«"^)+w(^"«)+w(-^^-^))/4. The 
f2("w) numerical flux given in Eq. (Eq) has now all 
the desired formal upwind properties both along the 
X and y directions, and reduces correctly to the 1-D 
limiting cases Eq. (|2^ ) and Eq. (|24|). The composition 
rule used here cannot be interpreted as a simple arith- 
metic average of independent 1-D Riemann solvers. In 
fact, for a discontinuity front with arbitrary slope an- 
gle around P, the 1-D flux formula, say Eq. (|2^) , still 
applies and gives the upwind contribution along the 
X characteristic modes. Near the y — yk+1/2 point, 
^x {b = L, R) is now a two-state flux function and 
upwinding has to be completed by taking into account 
also the — Aj, characteristic modes along the orthogo- 
nal y direction. By applying then 1-D flux upwinding 
as in Eq. ( p4| ) to the flux f7(w) = ^^(w) and by dis- 
carding quadratic terms, one recovers Eq. (EQ). This 
composition procedure taken in reverse order, start- 
ing now from Eq. (pj) , yields an identical result under 
the essential assumption of linearization. 

Finally, it is worth noticing that the viscous (resis- 
tive) model equation for the numerical flux function 
ft, consistent with Eq. (Eq), has the form 

showing how the dissipative term generalizes the clas- 
sical rjj = 77V X B term in Ohm's law of resistive 
plasmas. 

Having now completed the construction of the 
ri(w) upwind flux for the induction equations, it is 
possible to design the overall numerical procedure to 
integrate the MHD system. We summarize here the 
main computational steps: 

1. At each stage of the Runge-Kutta cycle, for 
given [u,A](t) data at time t, the averaged 
\Bx,By] staggered fields are evaluated first by 
Eqs. (0). 

2. Using then the {u, Bx,By)'^ data, all variables 
{u,Bx,By)'^ needed to compute the fiuxes de- 
fined in Eqs. (nsh are reconstructed at each cell 



boundary point (a;j+i/2j 2/fe) and conservative x- 
derivatives of the f flux can then be evaluated. 

3. The complementary procedure, now to recon- 
struct all variables {u,Bx,By)^ for the fluxes 
defined in Eqs. (M) at each {xj, yk+1/2) point, 
gives the conservative y-derivatives of the g flux. 

4. A final reconstruction to the {xj+1/2, yk+1/2) 
corner point is needed to compute the 17 (w) nu- 
merical flux in Eq. (|2^). This allows to integrate 
in time the vector magnetic potential A{t) by 
Eq. ^. 

The explicit representation of a divergence-free 
magnetic field via a vector potential has, among oth- 
ers, the advantage of an easy extension to three- 
dimensional configurations. In fact, in this case one 
has to discretize the constitutive relations 

B = VxA, B = [Bx,By,B,f 



in weak form to generalize Eqs. (M). The Ai, i = 
x,y,z components and the corresponding flux func- 
tions rii = [vx B]i arc now located at the 3-D cell edge 
points Pi, each centered along the corresponding i*'' 
coordinate but staggered with respect to the remain- 
ing two directions. The evolution equation Eq. (21) 
readily generalizes to 



d[A{t)] 



Pi 



dt 



f^-i(wp, 



and the consistency arguments introduced in 2D case 
can be applied to define the composition rules for each 
Cli — Qi{'w) scalar function since only two upwind 
directions, in turn, are now involved {i.e. each fli 
involves a 2D Riemann solver). 

3.3. Comments and discussion 

Some remarks are due in order to underline dif- 
ferences and analogies with other proposed MHD 
schemes, in particular with those presented by Dai & 
Woodward (1998) (DW), by Ryu et al. (1998) (RY) 
and by Balsara & Spicer (1999) (BS), all claiming to 
have "divergence-free preserving properties". These 
works are based on second order either Godunov or 
Roe-type schemes and a staggered discretization for 
the induction equation of the form of our Eqs. (EQ) is 
used. 

At a second order level the averaged variables 
[Bx]j+i/2,k and [By]jk+i/2 can be interpreted as 
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point valued variables at the same cell boundary 
points (which we label hereafter as {bx,by) fields, to 
conform to the DW and RY notation). In Sect. 3 we 
have shown that these staggered, divergence- free vari- 
ables, when used in momentum and energy equations 
avoid the effects of numerical monopoles. On the con- 
trary, in all the cited works, while evolving in time 
the staggered {bx,by) magnetic field, still interpola- 
tion and upwinding procedures are based on cell cen- 
tered {Bx, By)j^k variables and unwanted compressive 
V • B terms then necessarily set up. Moreover, this 
"duality" in the magnetic field representation is con- 
sidered to be unavoidable in the Godunov-type Rie- 
mann solvers formalism. For that reason, in the DW 
and BS papers in particular, much attention has been 
devoted to compare the different results produced by 
schemes advancing in time only the {bx,by) staggered 
variables, where (Bx,By) work as interpolated vari- 
ables, and schemes where also the {B^^By) cell cen- 
tered components are (independently) evolved (as for 
standard Godunov procedures for Euler equations). 
Conclusions are mainly drawn at a qualitative level, 
and only in the DW paper some numerical results on 
the V • B variable constructed with the (i?x. By) data 
are presented, showing the onset of significant, even 
of 0(1) size, residuals. 

In the present approach, we have demonstrated 
by analytical arguments that the compressive com- 
ponents arise either in cases where cell centered fields 
are evolved in time or they are simply given by inter- 
polation. In fact, at the leading second order inter- 
polation used in the cited papers, 

[Bx\j,k = ■^{[bx\j + l/2,k + [bx\j-l/2,k)^ 



field variables, depending on differentiation coordi- 
nates. The first set, given by the \bx{x^ ), 6y(, y)] vari- 
ables, which are continuous in the indicated coordi- 
nates, satisfy the divergence-free condition and does 
not have a characteristic representation. The second 
one, given by [Bx{,y),By{x,)], enters the character- 
istic space and can then have (transverse) disconti- 
nuities. The former quantities are advanced in time 
as staggered field data, while the latter can be recon- 
structed by interpolation either at a cell boundary or 
at a cell corner point. 

As a last point, we comment here on the way the 
flux for the induction equations is derived. In the 
DW and BS schemes, the Vt fluxes (or electric fields 
components) are constructed by a simple arithmetic 
averages in space and for each I-D upwind flux. This 
approach, which seems reasonable for second order ac- 
curacy, is not consistent since it does not reduce to the 
original I-D fluxes when the discontinuity front prop- 
agates along one of the coordinate axes. This draw- 
back led BS to introduce a rather empirical switch to 
select the dominant direction of front propagation. 
On the other hand, RY derives a formally correct 
flux, of the same form of our Eq. (|2g), by splitting 
the flux Q ^ VxBy — VyBx into two independent com- 
ponents {vxBy and VyBx), and 1-D independent up- 
windings along the x and y direction, respectively, 
have then been applied. However, this computational 
trick has no physical support, since the two charac- 
teristics spaces, spanned by the A^, and Ay, in the 
original MHD equations are both based on the com- 
plete J7 flux. 

4. Implementetion of a third order LF-CENO 
scheme 



[By]j,k - ■^{[by]j^k+l/2 + [by]3.k-l/2), 

one gets a 0{h^) (for smooth fields) residual when the 
V • B variable is evaluated by centered first deriva- 
tives. This compressive component associated to the 
{Bx,By) fields cannot be considered to be small (of 
the same order of the truncation error) , since it is easy 
to show that even for higher order (r > 2) interpo- 
lation the leading V • B = 0{h^) term never cancels 
out. 

Our analysis shows that cell centered fields are 
not really needed in upwind differentiation, even for 
Godunov-type schemes. In fact, the MHD system 
structure relies on two different kinds of magnetic 



We consider the 2-D MHD system in cartesian 
{x,y) coordinates, in the conservation form given by 
Eqs. ([l7|), to integrate the density, momenta and en- 
ergy variables, and by Eq. (g|) to integrate the vec- 
tor potential. The updated, line-averaged, magnetic 
fields {Bx, By) are defined, at each time step, by the 
geometrical relations Eqs. (|4|). 

We specify the general procedure outlined in the 
previous Sect. 3 by choosing high order reconstruction 
algorithms based on convex-ENO (CENO) method, a 
local Lax-Friedrichs (LLF) flux splitting for upwind- 
ing and a time integration step using a third order 
TVD Runge-Kutta scheme. 

All the indicated numerical ingredients are well 



12 



documented and tested in problems described by the 
Euler system of gas-dynamics (Shu & Osher 1988, for 
time integration, Shu 1997, for general ENO recon- 
struction and Liu & Osher 1998, for CENO method). 
It is then sufRcient to detail here the specific proce- 
dures allowing to extend CENO schemes to the MHD 
system and having relevance to the divergence-free 
properties of the magnetic fields. 

1. Among high order reconstruction algorithms, 
the recently proposed CENO method has the 
main computational advantage to avoid the time 
consuming characteristic decomposition of state 
variables, which is usually adopted in upwind 
schemes. 

To achieve this property, one consider first a 
TVD (monotone) second order accurate inter- 
polant for the scalar variable w{x) with data 
{wj}, where w denotes any component of the 



state vector w. In the 



Cj-l/2 



< X < X 



i-Hi/2 



range, three-point linear polynomials have the 
form 



L 



i'^)(^^- 



1 



J (x) =u;j + — [Aw] 



;+fe(x-Xj), k = 0,l 



where [Aw]i 



fj+i 



Wi 



In classical TVD 



schemes a unique intcrpolant Lj{x) = Wj + 
■^\P^^^w]j{x — Xj) with slope [P^^^wjj is cho- 
sen using min-mod (mm) limiters to assure 
nonoscillatory properties 



[D 



(i),„l . = 



w\i — mm 



{[Aw]j^i,[At 



The mm{a, b) function is defined, as usual, by 
mm{a, b) — sign{a)min[\a\, \b\], ah > 

and mm(a, b) — otherwise. 

Using the selected Lj{w) polynomial, the inter- 
polated values at the cell boundaries are then 
given by [w;(^)]j+i/2 = Zj^x^+i/a), [u;(^']j+i/2 = 

At any point of discontinuity as well at any 
smooth extrema of w{x) where the first differ- 
ences change sign, the TVD polynomial reduces 



to the constant state Lj{x) 



Clipping to 



first order accuracy at a function jump is un- 
avoidable in any polynomial based reconstruc- 
tion, while higher (r > 2) accuracy in smooth 
ranges can always be achieved by enforced TVD 



or ENO procedures. This improvement, how- 
ever, usually requires a preliminary decompo- 
sition of the Aw differences into characteristic 
modes, locally at each grid point. In the CENO 
method a higher order interpolation is main- 
tained only at smooth regions while first order 
polynomials are used at the function jumps, nei- 
ther cases requiring a characteristic decomposi- 
tion. 

For third order reconstruction, in particular, 
one has at disposal three quadratic intcrpolants 
in the Xj^i/2 ^ x < Xj^i/2 range 



Qf\x) ^w, + -{[Aw], + [Aw],^i) 



{X - Xi) 



where [A'^^^J^ = [A^ — Ai_i] and centering refers 
to the index i = j + k for k — —1,0, 1. The 
CENO selection procedure allows to construct 
the unique nonoscillatory interpolant 



,ix) 



i),„l.l£L^2l 



^. + -,[v^'W , 



which is closest to the TVD lower order inter- 
polant Lj{w). This is obtained by computing 
the three differences 

df^Qf-Lj, fc = -1,0,1 

at the X = a;j-i/2 or x = Xj+1/2 boundary 
point. In a smooth range all these distance in- 
dicators have the same sign and one can select 
then 

Only at a discontinuity point at least one in- 
dicator changes sign and in this case one takes 
Qj = Lj = Wj since V^^^ = V'--^^ = 0, chpping 
to a first order interpolation. 

2. We apply then this reconstruction procedure, 
first to each Uj,k scalar component of the u state 
vector, to recover the point values 



,(i)l 



+1/2, fe = Qj[u]{xj+i/2,yk) 
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needed to compute the {(\Vx,Bx) flux , and 
in a similar way to recover the Uj^k+1/2 values 
needed to compute the g(wy, By) flux. 

3. The definition of (B,j;,By) point values intro- 
duced in Sect. 3 requires the specification of the 
nonoscillatory V^^'' second derivative. In the 
CENO framework one simply takes 



V 



(2) _ 



(2) 



(2) a(2) 1 



2?(2) 
V 



.(2) A (2) A (2) 



We notice that this procedure returns the smoothest 
among the indicated three second numerical 
derivatives if the stencil of the involved first 
differences [Aj_2, ...., Aj+i] is monotone, while 
'D'2) = if the first derivative has a jump or a 
smooth extremum. 



For given [B^ 



j + l/2,fc: 



the reconstruction of the 



related divergence-free [i?x]j+i/2,fe point values 
is given by the implicit definition in Eqs. (p^ 



B,-7iI?i'^B,],+i/2,fe = [S,],+i/2,fc, 71 = 1/24, 



which we solve by using an explicit iteration 
procedure. By setting Bx ~ B^., at each 
{xj-^-l/2,yk) point, the sequence 



i?(")=B.+7iI?f[i?i"-^'] 



1,2,. 

-,(2) 



(27) 



is clearly rapidly convergent since Dx' is at 
most an 0{hx) quantity. In a similar way, we 
compute [i?j,]j,fc+i/2 by iterating 



bM 



By+j,Vi'^[Bi"-% 



5(") 



By. 



(28) 

In practical computations, for all the test prob- 
lems presented in the next section, we found 
that n = 5 is sufficient to assure V • B = to 
within machine accuracy both in the maximum 
and in the Li norm. 

The compute d \ Bx]i+i/2.k point values com- 
puted in Eq. ( P7|) enter now the {(wxyBx) flux 
as they stand, while the [By]j_^^i^2,k component 
of the Wx state vector needs a further inter- 
polation. Using then [i?y]j,fc+i/2 derived by 
Eq. (28), the cell centered [By]j^k field is first 
reconstructed by taking 



[By]jM 



2 ([-^2/] J, fe- 1/2 + [By]jM+l/2), 



\B.. 



y\j,k+l/2 



[B: 



y\],k+l/2 



;H'^B 



y\],k+l/2, 



to be finally interpolated at the (2^^-1-1/2, Vk) cell 
boundary point like the other Uj_k components 
of the Wx state vector. It is worth noticing 
that the cell centered [-By]j,fc values have not 
divergence- free properties since the {{wx,Bx) 
differentiation involves only the x coordinate. 

By symmetric arguments, the [i3j,]j.fe+i/2 field 
enters the g{wy,By) fiux as it stands, while 
[Bx]j+i/2,k has now to be interpolated using the 
cell centered values 



[Bx]j.k — -^{[Bx]j-l/2,k 



[By]j + l/2,k), 



\B._ 



x]j + l/2,k — [Bx]j+l/2,k - o [^2; Bx]j+i/2,k- 

4. The interpolated [vfx]j+i/2.,k or [wy]j^k+i/2 are 
represented as two-point left-right values along 
the relevant x or y coordinate, and an approx- 
imate Riemann solver has then to be specified 
to compute upwind fluxes. We have chosen the 
simple LLP flux composition, defined by 



f (w„ B,) = -[f (w^S B,) + f (wf s S,) 



--a4w,)(u«'-u^0, (29) 

and 



g(Wy,Ba 



[g{W^\By)+g{wf,By)] 



-ay{Wy){U''^~U-'^). 



(30) 



The scalar variable axi^x) is given at each 
{xj+i/2,yk) point by the largest of the A^, ma- 
trix eigenvalues As(wa;), and w^, is the arith- 
metic average of the w^; left-right states. In 
practice ax = \vx\ + Cf^, where c/^ is the fast 
wave speed along the x direction. Correspond- 
ingly, ay{\Vy) = Ivyl + Cf^ gives the largest 
eigenvalue of the Ay matrix based on the arith- 
metic average of the left-right states of w^ at 
the {xj,yk+i/2) point. 



5. Differences of flux values given in Eqs. (E!*, (30) 



provide only second order accurate derivatives, 
even if the reconstructed flux arguments share 
a higher accuracy order. To keep third order 
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in derivative approximations, we construct the 
primitives 



rj+l/2,fe 



= [f-7i2?i')f] 



j + l/2,fcj 



k],k+l/2 = [g - llT^y 's]j.k+l/2 

thus completing the integration scheme of Eqs. (Il 
for the six-component state vector u^ j.. 

6. The 0(w) flux variable needs a proper upwind- 
ing procedure, as shown in Eq. (Pq). In a LLP 
scheme one has 

fl{wp) = 

- ■^[ax{^)SxBy - ay{w)SyB^]p, (31) 

where all the arguments {vx, Vy, Bx, By) are first 
interpolated at a common P — {xj+1/2, yk+1/2) 
point. 

7. Finally, for time integration, a three-step Runge- 
Kutta algorithm provides the overall third order 
accuracy of the LF-CENO scheme. 

5. Numerical results 

The proposed numerical problems are mainly con- 
cerned with the divergence-free property, which on 
numerical side entails two main aspects: 

1. the existence of a vector potential A{xj^i/2-, yfc+1/21 1) 
as a continuous function at all times, assuring 
regular field lines topology (only corners are al- 
lowed) , 

2. a vanishing Db = [V • 'E]j^k where field deriva- 
tives are computed using the same field compo- 
nents and the same difference algorithms as in 
dynamical flux calculations. 

Since this form of validation has no counterpart in 
other proposed numerical works, comparisons with 
published data will cover necessarily rather qualita- 
tive aspects. Beside the divergence-free condition, the 
numerical results give also indications on the resolu- 
tion properties, as well as on the stability and relia- 
bility of the MHD code described in the Sect. 4. 

Finally, we remark that in our code the A{x, y) vec- 
tor potential refers only to the nonuniform {Bx,By) 



fields, since constant initial components are trivially 
preserved in time. Therefore, for problems having 
constant components {Bq^^ B^y), the evolved A(t) 
field is now defined by 

Bx = Bf)x + dyA, By — Boy — dxA, 

replacing the original Eqs. (Eoh of Sect. 2. 

5.1. Shock-tube tests 

We first consider 1-D Riemann problems (using a 
full 2-D grid) to check for resolution properties of the 
proposed scheme. To that purpose, it is necessary 
to take into account that high order schemes are not 
well suited for shock-tube problems where lower order 
characteristic-based schemes are optimal, instead. 

We consider three problems documented by Ryu 
and Jones (1995) in their Fig la. Fig 2a and Fig 5a, 
which we here label RJl, RJ2 and RJ3, correspond- 
ingly. In all the indicated cases a uniform grid with 
Nx = 400 grid points, a grid size L^; = 1 , an adiabatic 
index 7 = 5/3 and a CFL number c = 0.8, are used. 
In all numerical tests presented here, the parameters 
7 and c will always retain the same values. 

In RJl, the initial conditions for the state vector 
w{x) = [p,Vx,Vy,Vz,By,Bz,p\'^ are defined by 

w^ = [1, 10, 0, 0, 5Bo, 0, 20]^, 

w^ = [1, -10, 0, 0, 5^0, 0, 1]^, 

and by a constant Bqx — 5i?o- Here left states refer 
to 2: < 0.5 and right states to x > 0.5. The unit 
magnetic field is Bq = 1/^/4tt. In Fig. || the evolved 
variables {p,p, By,Vx) are shown at time t = 0.08, as 
in the referenced RJ paper. 

In the RJ2 test initial conditions are defined by 
w^ = [1.08, 1.2, 0.01, 0.5, 3.6B0, 2Bo, 0.95]^, 

w^ = [1, 0, 0, 0, 4Bo, 2Bo, if, 

and the constant magnetic field is now Bqx = 2Bq. 
The evolved variables (p,p. By, Vy, Bz,Vz) at time t = 
0.2 are shown in Fig. 0. 

Finally, the RJ3 problem, with initial data 
w^ = [1,0,0,0,1,0,1]^, 

w^ = [0.125,0,0,0,-1,0,0.1]^, 

and Bqx = 0.75, is illustrated in the Fig. ^ for t = 0.1. 
This is already considered a classical test, related to 
the presence of a compound wave (Brio & Wu 1988). 
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As can be seen, the plotted results reproduce well 
all the main expected features and compare with the 
corresponding results obtained with higher grid res- 
olutions and more elaborate Riemann solvers (Ryu 
& Jones 1995, Jiang & Wu 1999). Postshock os- 
cillations, which are always produced in any shock- 
capturing scheme (Arora & Roe 1996), appear here 
with vanishing amplitudes behind the fast moving 
shocks but have significant sizes near the slow shocks 
and near the expansion wave on the right hand side 
of Fig. ^. At present, to our knowledge, no general 
cure has been envisaged to suppress entirely this un- 
physical wave noise, which can then only be reduced 
by adding numerical viscosity. In this sense, the ob- 
served oscillations allow to estimate the implicit nu- 
merical viscosity of our CENO-LF scheme to be some- 
how intermediate between the lower order TVD code 
of Ryu & Jones (1995) and the Weno-LLF MHD code 
of Jiang & Wu (1999). Other hmiters have also been 
tested (van Leer, superhee, and so on), with no signif- 
icant improvements. In any case, we want to stress 
again the point that high order schemes not based on 
characteristics decomposition, like our code, are not 
particularly designed to handle Riemann problems, 
where a lower order scheme may be a better choice. 

To test the code for 2-D cases, we have run the 
previous Riemann problems RJl and RJ2 with struc- 
tures propagating along the main diagonal of a com- 
putational box with sizes Lx = cos a, Ly = sin a, 
where a = 7r/4. In this way the diagonal has a unit 
size L = 1 and h = 1/Nx is the size of the cell diago- 
nals. Initial conditions are then assigned to the state 



vector w(^) 
ordinate ^ = 



[p,V^,V^,V;, 



Brf,Bz 



p] along the co- 



xcosa + ysina, with now B^ = B{ 



being constant. Boundary conditions are specified by 
imposing the continuity of all variables along the tra- 
verse direction 77 — ycosa — xsina, extended to the 
X < Q, X > Lx and to the y < 0, y > Ly sides. We 
used Nx — Ny — 256 grid points and we found that 
this grid spacing is hardly sufficient to recover the 
main flow structures. 

The evolved variables w(^) are shown in Fig. ^ 
for the rotated 1-D Riemann problem of Fig. |^, and 
in Fig. H for the rotated 1-D Riemann problem of 
Fig. 0, at corresponding times. The plotted results 
compare to the ones presented by Ryu et al. (1998) for 
the same Riemann problems. The fact that the small 
oscillations observed in the corresponding 1-D cases 
are now less apparent, is due to the higher numerical 
dissipation produced by the lower resolution. 



In a 2-D shock-tube problem, the divergence-free 
condition can be simply expressed by a constant B^ 
field, i.e. by the dB^/d^ = relation along both the 
^ and r] coordinates. However, if i?^ is constructed 
using the cell centered [Bx,By]jk fields, this conser- 
vation law is poorly verified, as can be seen in Fig. pk, 
where the numerical derivative A^B^/h for the RJ2 
test is plotted. On the other hand, using the vector 
potential, point values of the B^ field are properly 

defined by 

dA 

and the divergence-free condition results if A (as well 
as the other dynamical variables) do not depend on 
the T] coordinate. This is documented by the 2-D 
structure shown in Fig. gb. The corresponding nu- 
merical derivative 



1 A^Ar,A 



has now a maximum size of ~ 10 ^ . 

5.2. Slow v^rave steepening and shock forma- 
tion 

Nonlinear wave steepening from continuous initial 
data is a main feature of compressible flows. In the 
MHD case this problem has also interesting astro- 
physical aspects for the study of intermediate shocks 
(shocks coupled to expansion waves), already encoun- 
tered in the previous shock-tube test RJ3. On the 
computational side, wave steepening is significant for 
high order schemes, whose small numerical dissipa- 
tion may model weakly resistive plasmas. Reference 
results are given by Wu (1987) for physical setting in 
resistive MHD and by Dai & Woodward (1998) and 
Jiang & Wu (1999) for numerical testing. 

Here we consider initial data defined by a (smooth) 
slow wave front propagating along the transverse ^ 
axis, with a slope angle a = -k /Q with respect to the x 
axis. The initial conditions are defined by the charac- 
teristic differential equations (the prime here denoting 
a ^ derivative) 



P = 



a2 - c2) 



, P 



2 I 
a p , 



9? 



—p 
p 






relating {p, p, g^, g^) to the _B^ field along the ^ co- 
ordinate. The variables a, c^, c/ denote the sound, 
slow and fast wave speeds, respectively. We choose as 
initial profile i3j;(0 = sin(27r^) and a cartesian box 
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with Lx — 1/cosa and Ly — 1/sina, so that peri- 
odic boundary conditions can be applied along the x 
and y cartesian coordinates. 

In Fig. the ^ profiles of the variables (p, p, i?,,, v^ ) 
are shown at time t = 1, when a shock train is al- 
ready formed. The corresponding 2-D plots of the 
pressure p{x,y) and of the vector potential A{x,y) 
are also shown in the Fig. K to check for accurate rj 
independence of the flow variables. As for the pre- 
vious shock-tube test of Fig. 0, a vanishing numer- 
ical V • B comes from the drjA ~ condition. In 
Fig. Ha a surface plot of the variable Db is also shown, 
giving a value \DB\max — 10^^ for the residual nu- 
merical monopoles, while the B^ component based 
on {B^, By)jk values shows a much higher derivative 
0(10^^), as can be seen in the Fig. ||b. 

5.3. The Orszag-Tang MHD vortex problem 

A well known model problem to study the transi- 
tion to MHD turbulence is provided by the so-called 
Orszag-Tang vortex, which has been extensively stud- 
ied in its compressible version (for low Mach numbers) 
by many authors using spectral methods. For initial 
Mach numbers M >\ this is also a valuable test for 
upwind codes, and it has been used for almost all the 
latest schemes (Zachary et al. 1994, Dai & Woodward 
1998, Ryu et al. 1998, Jiang & Wu 1999). The ref- 
erenced Orszag-Tang system is defined by the initial 
conditions 



Bx/Bq = -sin27r2/. 



Vy = sin 2-1^ X, 



By = i?o sin47ra;. 



p ^ {(i/2)Bl, p = 7p, 

where Bq — l/\/4n and [3 — 2^ for the usual 7 = 5/3 
value. The initial flow is the given by a velocity vortex 
superimposed to a magnetic vortex, with a common 
(singular) X-point, but with a different modal struc- 
ture. This configuration is strongly unstable, giving 
rise to a wide spectrum of propagating MHD modes 
and shock waves (here the initial Mach number is 
A/ = 1), and to the transformation of the initial X- 
point to a current-sheet triggering the reconnection 
process. 

For this test we have chosen a unit grid L^ = Ly = 
1 with Nx — Ny = 192 collocation points. In Fig. |l^ 
we present the pressure p and potential A isocon- 
tours at t=0.5, showing the good qualitative agree- 
ment with the other published works. In particular, 
1-D profiles of the p{x) variable at y = 0.4277 (up- 
per plot) and at y — 0.3125 (lower plot) are shown 



in Fig. |T^ for a more detailed comparison with the 
corresponding plots given, respectively, by Ryu et al. 
(1998) and Jiang & Wu (1999). The latter reference 
also contains some quantitative estimate of how sig- 
nificant magnetic monopoles may affect the computed 
solutions, thus producing numerical instabilities in a 
long time evolution. 

The monopole distribution Db (not plotted here) 
show only a few enhanced values ~ 10~^, thus as- 
suring vanishing V • B condition, in the long time 
computations we found no evidence of negative pres- 
sure nor other unphysical behaviors. To check this 
point in more detail, we plot the magnetic field lines 
at t = 3 in Fig. 12, showing how the regularity of 
fieldlines is well preserved in time. This essential fea- 
ture of the computed divergence-free magnetic field 
allows to reproduce typical magnetic phenomena, like 
the topology change induced by a vanishing resistiv- 
ity (here modeled by a low numerical diffusivity). In 
fact, by comparing the latter A distribution with the 
former shown in Fig. 10, it is apparent how the ini- 



tial magnetic islands around the X-point merge by 
reconnection. 

5.4. Strong blast wave in free space 

The following two numerical tests concern the for- 
mation and propagation of strong MHD discontinu- 
ities in a 2-D domain. These model problems are rep- 
resentative of many astrophysical phenomena where 
the magnetic energy has relevant dynamical effects. 
In numerical schemes having poor divergence-free 
properties, the (possible) onset of spurious solutions 
and of a negative gas pressure is clearly enhanced in 
these physical regimes, since the magnitude of numer- 
ical monopoles increases along with the background 
magnetic pressure. This problem has been discussed, 
in particular, by Balsara & Spicer (1999), where also 
some quantitative estimate of the numerical V • B pro- 
duced by Godunov schemes has been documented. 

The first test problem concerns the explosion of a 
circular dense cloud in a magnetized, initially static 
region. Here we take again a square domain with 
Nx = Ny — 192 grid points. Initial conditions are 
specified by filling a circular region located at the 
center and radius tq — 0.125 with a hot gas having 
p = 100. The background static fiuid is characterized 
by P = P = 1 and Bqx — 10. 

In Fig. |l^ we show the density p, the magnetic po- 
tential A^ the magnetic pressure (i?^ -|-i?^)/2 and the 
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solenoidal variable Db distributions at time t ~ 0.02, 
which is already representative of the generated com- 
plex flow structure. In particular, the plotted re- 
sults show the well preserved initial axial symmetries 
(around both the y = 0.5 and the x = 0.5 axis) as well 
as the regularity of magnetic fieldlines. As we can see, 
the resulting numerical Db variable has an isolated 
peak with magnitude 10"'^ and otherwise vanishing 
sizes < 10^^. 

5.5. The fast rotor problem 

In the Balsara & Spicer (1999) work, a model prob- 
lem to study the onset and propagation of strong 
torsional Alfven waves, relevant for star formation, 
has been presented and analyzed. Following this ref- 
erence, we have runned the same problem using a 
square unit computational box and Nx = Ny — 240 
grid points. For propagating structures not intersect- 
ing the box boundaries, periodic conditions can be 
applied. Initial conditions are specified by a rapidly 
rotating cylinder (the rotor) with center at the x = 
0.5, y — 0.5 point and radius r = 0.1. The rotor 
has (initial) density p = 10, angular velocity lu — 20 
and it is embedded in a static and uniform fluid with 
p = p = 1 and Bqx = 2.5/-\/7r. 

The flow pattern evolved at time t = 0.18 (just be- 
fore the shocked flow reaches the boundaries) is shown 
in Fig. |lj, representing as in Fig. ^ the space dis- 
tribution of the density p, the magnetic potential A, 
the magnetic pressure {B'^ + By)/2 and the solenoidal 
variable Db, that keeps everywhere below 4 x 10~^. 
Even if a lower resolution than in the referenced paper 
is here adopted, and no smoothing has been applied 
to the initial density and rotation velocity discontinu- 
ities, the numerical results give convincing evidence 
on how a higher order scheme provides accurate and 
well resolved profiles even when strong discontinuities 
develop. 

6. Conclusions 

We have introduced a general method to adapt 
upwind schemes developed for Euler system to the 
corresponding MHD system in order to assure the 
divergence-free condition. The proposed approach 
can be applied to existing MHD codes as well as to 
any higher order extensions. 

The use of a staggered collocation for the magnetic 
field components entering the V • B variable and of 
the related magnetic potential A are well known gen- 



eral premises to represent a numerical divergence-free 
magnetic field at the second order accuracy level. We 
have thus introduced proper algorithms to extend this 
representation to higher orders and to formulate up- 
wind flux derivatives using only the divergence-free 
variables, in order to avoid the onset of numerical 
monopoles in the momentum and energy equations. 
Moreover, by taking into account consistency argu- 
ments, we have proposed a new formulation for the 
upwind flux for the induction equations. 

As an application, we have constructed a sim- 
ple and efllcient third order LF-CENO based MHD 
code running in multidimensional systems. This code 
appear well suited for many astrophysical problems 
where, beside strong shocks, reconnection phenom- 
ena, complex wave patterns and turbulence develop, 
as confirmed by the several numerical tests here pre- 
sented. 

The authors would like to thank Marco Velli for 
many helpful discussions and for his support in com- 
pleting this work. 
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Fig. 1. — The indicated variables at time t — 0.08 
for the 1-D Riemann problem RJl, using N^ = 400 
grid points. 

Fig. 2. — The indicated variables at time t = 0.2 for 
the 1-D Riemann problem RJ2, using N^ — 400 grid 
points. 

Fig. 3. — The indicated variables at time t = 0.1 for 
the 1-D Riemann problem RJ3, using N^ = 400 grid 
points. 

Fig. 4. — The same shock-tube problem as in Fig. |^, 
now along the main diagonal ^ of a 2-D square com- 
putational box with (256 x 256) grid points. 

Fig. 5. — The same shock-tube problem as in Fig. 0, 
now along the main diagonal ^ of a 2-D square com- 
putational box with (256 x 256) grid points. 

Fig. 6. — Left panel (a): the ^ derivative of the B^ 
field in the 2-D shock-tube problem of Fig. Q Right 
panel (b): the corresponding isocontours of the mag- 
netic potential A, which clearly show the 77-invariance. 



t = 3. The reconnected central magnetic island is 
clearly shown. 

Fig. 13. — The indicated variables at time t = 0.02 
for the blast wave problem. A unit computational 
box is used with (192 x 192) grid points. 

Fig. 14. — The indicated variables at time i = 0.18 
in the rotor problem. A unit computational box is 
used with (240 x 240) grid points. 



Fig. 7. — The indicated variables for the 2-D slow 
wave problem, along the ^ coordinate and at time t = 
1. The slope angle is a = tt/6 and the computational 
box has (192 x 192) grid points. 

Fig. 8. — The pressure (left) and magnetic potential 
(right) isocontours for the slow wave problem of Fig. ^ 
at the same time t = 1. 

Fig. 9. — Left panel (a): A surface plot of V • B for 
the slow wave problem at time t = I. Right panel (b) : 
The numerical derivative of the B^ field, computed 
with cell centered {Bx,By) data. 

Fig. 10. — The pressure (left) and magnetic potential 
A distribution (right) in the Orszag-Tang problem at 
time t = 0.5. The unit square computational box has 
(192 X 192) grid points. 

Fig. 11. — The 1-D pressure distribution for the same 
problem as in Fig. |l^ along a cut at y = 0.4277 (upper 
plot) and at y = 0.3125 (lower plot, where a proper 
normalized pressure is shown, to compare with the 
Jiang & Wu (1999) data). 

Fig. 12. — The magnetic potential A distribution 
for the same problem as in Fig. fo at a later time 



20 



Q. 



5 r 

4 

3 h 
2 r 

1 



0.0 0.2 0.4 0.6 0.8 1.0 




m 




> 



1b 


: 


: 


10 


I O 


-_ 


5 


: o 


"■_ 









-5 


- 


- 


10 


- 




1,S 







0.0 0.2 0.4 0. 



0.8 1.0 



21 




2.U 




' 




- 


1.8 


- 


r 

o 

J 


i. 


- 
















o 






"^ _ 


,z 


o 






: 


1.0 


» 






4_ 






n R 




, 


1 


- 



U.J> 






' 




E 


n.? 


L 




* 




J 








1* 




- 








i« 




= 




















1 




- 






^ 


Xm 


•^^ 


- 












0m 


I) (1 
















-- 




























^ 
























If 


I 


9 










-- 



0.1 

m" o.e 

o.e 



I o I / ^ J 

^ O O Oh 

t \ / j 



0.2 
0.0 

-0.2 




22 



1.2 
1.0 

0.8 

Q. 0.6 

0.4 

0.2 

0.0 





0.2 




U 



%• 



0.4 0.6 0. 

X 



1.2 
1.0 

0.8 

0.6 
0.4 
0.2 

0.0 

-0.2 






0.2 0.4 0.6 0. 



1.0 



CD 




> 



0.0 0.2 0.4 0.6 O.l 




23 



Q. 



5 r^ 

4 

3 h 
2 r 

1 - 



0.0 0.2 0.4 0.6 



O.f 



1.0 




0.0 0.2 



m 




0.0 0.2 



0.4 



.6 0.8 



> 




24 



1.8 










- 


1.6 


- 




f 


\ 


- 


1.4 






1.2 


1 








o 1 












w 


1.0 


7 




OR 










- 



a.u 







, 


. 


1.8 


r 






-_ 




_/ ^^ 


o - 


1.4 






1 P 




" 




J 


1.0 


1 


, 




%_ 








R 




1 1 


, 1 


- 



0.8 L_ 

0.0 



0.3 




' 


' 


E 


0.2 






Sb 


J 










^ 


0.1 




°*— ^« 










j- 






o I 




r~~ 


I 


-0.2 


r 






— 


-0.3 








E 



0.90 
0.80 


h 








J 


o 




0.70 


f^ 




m 0.60 


L 




o 


o 

-J 


^ 


0.50 






0.40 


h 








i 


30 










E 





1.0 










: 




0.8 


- 




«H : :-4^ 




- 








o 


, 




- 




0.6 


- 








1 








>" 


0.4 


- 




V. 


^— ^ 


- 




0.2 


r 






\ 


^ 




0.0 

-n ? 








K— ™. 


S— 









0.4 0.6 



25 



Magnetic potentia 




0,4 0.6 



0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 



26 




0.0 



0.5 



1.0 



1.5 



2.0 





2.0 



0.05 


0.00 


-0.05 


r -0.10 


-0.15 


-0.20 


-0.25 







27 



Pressure 



Magnetic potentia 



>^ 1.0 




0.0 0.2 0.4 0.6 0.8 1.0 

X 



>^ 1.0 




0.0 0.2 0.4 0.6 0.8 1.0 

X 



28 





29 



Pressure T = 0.5 



Magnetic potential T=0.5 





0.0 0.2 0.4 0.6 0.8 1.0 



0.0 0.2 0.4 0.6 0.8 1.0 



30 




0.00 t 



0.0 



0.2 



0.4 



0.8 



.0 





4 


z 




z 




3 


< 

Mill 











s<«*e>o " ^^V 


^^^KXXX^X^^^ ^ 


CL 


2 


^ 







] 


E 


^ ^p,:^ 






= 


«^ 


^--^^^^^^'X^^^ = 




n 


z 




z 



0.0 



0.4 



0.8 



.0 



31 



aqnetic potential T = 3 



>^ 




32 



Density 



Magnetic potentia 



0.8 - 



0.6 - 



0.4 ■ 



0.2 - 



0.0 [ 




0.8 



0.6 



0.4 



0.2 



0.0 




0.0 0.2 0.4 0.6 0.8 1.0 



0.2 0.4 0.6 0.8 1.0 



Magnetic pressure 



0.8 - 



0.6 - 



0.4 



0.2 - 




0.0 0.2 0.4 0.6 O.i 




33 



Density 



Magnetic potentia 



0.8 - 



0.6 - 



0.4 ■ 



0.2 - 



0.0 [ 




0.8 



0.6 



0.4 



0.2 



0.0 



0.0 0.2 0.4 0.6 




Magnetic pressure 




0.0 0.2 0.4 0.6 0.8 1.0 




34 



